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Abstract  This  paper  provides  a  novel  incremental  multigrid  strategy  for  the  equations  of  fluid  dynamics. 
The  (time  dependent)  governing  equations  are  discretized  in  time  by  means  of  a  two  level  implicit  Euler 
scheme  and  linearized  using  Taylor  series  and  the  incremental  (delta)  form  of  Beam  and  Warming.  The 
coefficients  and  the  right  hand  side  of  the  resulting  linear  systems  are  evaluated  always  at  the  finest  grid 
level,  whereas  the  (delta)  unknowns  are  computed  (approximately,  by  a  single  relaxation  sweep)  on  a 
sequence  of  coarser  meshes.  At  every  grid  level  the  computed  deltas  are  interpolated  up  to  the  finest- grid 
level  and  used  to  update  the  solution,  as  well  as  the  coefficients  and  the  right  hand  side  of  the  linear 
systems.  The  process  is  repeated,  sweeping  all  grid  levels  successively,  until  a  satisfactory  convergence 
criterion  is  met.  The  validity  of  the  proposed  approach  is  demonstrated  by  solving  a  simple  linear 
problem  and  the  vorticity-stream  function  Navier-Stokes  equations,  using  line  relaxation  methods  as 
smoothers,  and  the  lambda-formulation  Euler  equations,  in  conjunction  with  a  simple  explicit  smoother. 
In  all  cases,  the  proposed  multigrid  strategy  provides  a  considerable  efficiency  gain  over  the  corresponding 
single-grid  methods. 
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Introduction 

Id  the  last  years,  tbe  author  has  been  working  on  developing  and  improving  implicit  numerical  methods 
for  solving  the  Navier-Stokes  equations  and  the  lambda-formulation  Euler  Equations,  mostly  for  the  case 
of  steady  flows.  The  basic  framework  of  all  numerical  schemes  developed  by  the  author  and  his 
coworkers  (see,  e.  g..  [1-7])  is  the  following:  the  time  dependent  governing  equations  are  discretized  in  time 
by  a  two  level  implicit  Euler  time  stepping  and  linearized  using  Taylor  series  and  the  incremental  (delta) 
approach  of  Beam  and  Warming  [8];  the  resulting  large  block-sparse  linear  system  obtained  at  every  time 
step  is  reduced  to  a  series  of  smaller  block-tridiagonal  systems  by  approximate  factorization  [1-5]  or 
simple  "mutilation*  [6,7];  these  systems  are  solved  very  efficiently  by  means  of  standard  block-tridiagonal 
Gaussian  elimination  [9];  the  solution  is  updated  and  the  process  is  repeated  until  a  satisfactory 
convergence  criterion  is  met. 

The  rationale  behind  such  a  general  approach  (see  also  [8,  10-12])  is  that  a  single  relaxation-like  procedure 
(requiring  to  solve  only  block-tridiagonal  systems)  is  used  to  relax  for  both  the  nonlinearity  of  the 
governing  equations  and  the  extra-tridiagonal  terms  in  the  linear  systems  arising  at  every  time  level, 
leading  to  a  globally  efficient  procedure.  A  major  benefit  of  employing  the  incremental  approach  [8]  is  in 
that  a  deferred  correction  strategy  [13,14]  can  be  implemented  very  elegantly  and  straightforwardly  (see, 
e.g.,  [4]):  two-point  first-order-accurate  upwind  differences  are  used  for  the  advection  terms  in  the  LHS 
(left  hand  side)  implicit  operator,  whereas  second-  or  higher-order-accurate  differences  are  used  for  the 
RHS  (right  hand  side)  steady-state  equations,  which  are  evaluated  at  the  old  time  level  explicitly.  In  this 
way,  the  linear  systems  to  be  solved  at  every  time  step  are  guaranteed  to  be  diagonally  dominant  so  that 
improved  convergence  rate  and  stability  are  obtained  without  lowering  the  accuracy  of  the  final  steady- 
state  solution,  in  general,  and  classical  point  and  line  relaxation  methods  have  regained  a  competitive 
edge  with  respect  to  more  sophisticated  approximate  factorization  methods  [6,7,16],  in  particular.  Is  is 
noteworthy  that  from  a  fluid  dynamics  point  of  view  (for  the  case  of  the  Navier-Stokes  equations),  the  use 
of  such  a  deferred  corrector  approach  is  a  means  of  stabilizing  tbe  pseudo-transient  by  adding  artificial 
viscosity,  which  is  proportional  to  a  time  derivative  and  thus  diminishes  (and  eventually  vanishes)  when 
approaching  (and  finally  reaching)  the  steady  state. 

The  aforementioned  "CFD  philosophy"  has  been  reasonably  successful,  especially  from  the  point  of  view 
of  its  general  applicability  to  the  compressible  Navier-  Stokes  [11,8,6]  and  Euler  [1,2,7]  equations  as  well 
as  to  the  incompressible  vorticlty-stream  function  (and  energy)  equations  [3-6]  governing  the  motion  of 
two-dimensional  viscous  flows  (including  buoyancy  effects). 

However,  the  convergence  rate  of  all  of  these  methods  Invariably  deteriorates  when  the  mesh  is  refined,  so 
that  for  very  complicated  problems  requiring  a  very  fine  computational  grid,  they  tend  to  lose  their 
■competitive  edge". 

The  aim  of  this  paper  is  to  remove  such  a  limitation  by  developing  a  simple  incremental  multigrid 
strategy  especially  "custom-tailored"  for  all  of  these  methods,  so  as  to  improve  their  convergence  rate 
especially  at  the  later  stages  of  the  relaxation  process. 

There  are  two  complementary  and  somewhat  antithetic  ways  of  looking  at  a  multigrid  strategy  as  a 
means  for  obtaining  a  desired  numerical  solution  efficiently  [16]:  the  first  way  is  to  consider  every  finer- 
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grid  computation  as  a  means  for  Improving  tbe  accuracy  of  a  solution  cheaply  obtained  on  a  coarser 
mesh;  the  second  way  is  to  regard  every  coarser-gTid  computation  as  a  means  for  reducing  a  given 
frequency  component  in  the  residual  as  effectively  as  possible.  Of  course,  every  articulate  multigrid 
strategy  combines  both  the  aforementioned  aspects.  Insofar  as  it  makes  use  of  smoothing  processes  at  all 
grid  levels  (except  maybe  at  the  coarsest  one)  as  well  as  of  prolongation  and  restriction  operators  among 
the  various  grid  levels.  However,  by  considering,  for  example,  three  very  successful  multigrid  methods 
employed  for  solving  flow  problems  [17-19],  it  is  clear  that  tbe  method  of  Ghia  et  al.  [17]  is  more  typical 
of  the  first  point  of  view.  Insofar  as  a  converged  solution  is  obtained  on  successively  finer  grids  as  a 
starting  point  for  smoothly  and  efficiently  achieving  the  sought  flnest-grid  solution;  the  methods  of  Ni 
[18]  and  Jameson  [19]  are  Instead  more  typical  of  the  second  point  of  view,  insofar  as  they  only  obtain  the 
flnest-grid  solution  and  use  the  coarser  meshes  Just  to  annihilate  the  residual  as  fast  as  possible. 

The  present  approach  has  been  devised  to  improve  the  convergence  rate  of  the  various  numerical  methods 
developed  and  employed  by  the  author  to  solve  fluid  dynamic  problems  [1-7].  From  the  very  nature  of 
these  methods,  it  is  clear  that  the  proposed  multigrid  approach  is  more  similar  to  those  of  Nl  [18]  and 
Jameson  [19].  As  a  matter  of  fact,  the  present  methodology  goes  a  step  further  towards  employing  a  single 
relaxation  sweep  on  each  grid  to  obtain  the  final  steady  state  solution  on  the  finest  grid  as  fast  as 
possible,  insofar  as  tbe  incremental  variables  computed  from  each  coarse-grid  relaxation  sweep  are 
immediately  used  to  update  the  flnest-grid  solution;  only  the  flnest-grid  solution  is  ever  computed, 
whereas  the  deltas  are  evaluated  sequentially  on  a  series  of  coarser  grids  by  means  of  a  single  relaxation 
sweep,  immediately  Interpolated  back,  up  to  the  flnest-grid  level  and,  finally,  summed  to  the  flnest-grid 
solution,  which  is  therefore  updated  after  every  relaxation  sweep. 

The  proposed  methodology  is  described  in  the  particular  case  of  the  vorticity-stream  function  Navier- 
Stokes  equations,  in  the  next  section,  and  Is  then  applied  to  three  different  sets  of  equations,  namely;  the 
Laplace  equation,  the  aforementioned  Navier-Stokes  equations  and,  finally,  the  lambda-formulation  Euler 
equations. 

Methodology 

The  proposed  multigrid  strategy  is  described  here  for  convenience  for  the  particular  case  of  the  vorticlty 
stream  function  Navier-Stokes  equations: 

wt  +  Wx  *  ^XWy  '  (WXX  +  Wyy)/Re  ~  0  U) 

^  +  V’y,  +  W  -  Vt  =  0  (2) 

In  eqns  (1-2)  Re  is  the  Reynolds  number,  -w  is  the  vorticlty,  \l>  is  the  stream  function,  x  and  y  are  the 
standard  Cartesian  coordinates,  t  is  the  time,  subscripts  indicate  partial  derivatives  and  a  relaxation-like 
time  derivative  has  been  added  to  the  stream  function  equation  to  parabolize  it.  Eqns  (l,  2)  are 
discretized  and  linearized  in  time  using  a  two  level  implicit  Euler  scheme  and  the  delta  approach  of  Beam 
and  Warming  [8]  to  give: 

A^/At  +  t/>yn  Au,"  +  -  wyM«/>xH  -  (Aw^  +  4wyyH)/Re  = 
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-  (^n  u\  +  (V-xDwn)y  +  (WjocD  +  u,yyD)/Re  =  RES(u,)D  (3) 

AijP/At  -  AipJ*  -  AiPyjH  -  AuH  =  ^  +  V>„n  +  wD  =  RES(^)“  (4) 

In  eqns  (3,  4)  At  is  the  time  step  (which  can  be  different  in  eqns  (3)  and  (4)  and  can  also  vary  at  every 
iteration).  Aw  —  wn+1  -  wn,  the  superscripts  n+1  and  n  Indicating  the  new  and  old  time  levels  tn+1  and 
tD.  etc.  and  the  superscript  H  Indicates  the  mesh  on  which  the  superscripted  variables  are  evaluated 
(notice  that  the  superscript  h.  indicating  the  finest  mesh  employed  in  the  computation,  is  always  omitted, 
for  convenience);  finally,  the  conservative  form  of  the  convection  terms  is  used  in  the  RHS  of  eqn  (3), 
since  it  has  been  shown  to  provide  improved  accuracy  over  the  standard  form,  especially  at  high  Re 
values  (see,  e.  g.,  [17,  8]). 

Eqns  (3*4)  are  discretized  in  space  using  second-order-accurate  central  differences  throughout  (except  for 
the  incremental  convective  terms  of  the  vorticity  equation  in  the  LHS  of  eqn  (3),  which  are  approximated 
using  first-order-accurate  upwind  differences),  so  that  a  large  linear  system  is  obtained,  characterized  by  a 
2x2  block-pentadiagonal  matrix.  Such  a  linear  system  is  then  solved  approximately  (by  a  single-sweep 
relaxation  method  requiring  to  solve  only  block-tridiagonal  systems).  The  process  is  started  on  the  finest 
grid,  with  A(J*  —  Aui,  Aip*1  —  Aip,  the  solution  {ipD,  wn)  is  updated  and  Eqns  (3,  4)  are  then  solved  on 

successively  coarser  grids  (H  —  2h,  4h . ),  until  the  flnest-grid  residual  is  reduced  to  a  suitable  small 

value.  In  more  detail,  at  every  grid  level  H.  the  following  steps  are  required  by  the  proposed  multigrid 
strategy:  a)  the  coefficients  of  eqns  (3,  4)  are  evaluated  from  the  flnest-grid  (updated)  solution  ipn,  wB, 
locally  at  the  H-mesh  gridpoints  (which  amounts  to  injecting  the  h-mesh  coefficients  into  the  H  mesh);  b) 
the  RHS  of  eqns  (3,  4)  are  replaced  by  the  corresponding  collected  residuals,  given  as: 

RESH(V>,  w)  =  ChH  RES(V>.  w)n  (5) 

where  Ch  indicates  a  suitable  collection  operator  from  the  finest  grid  h  to  the  current  grid  H;  c)  eqns  (3, 
4)  are  solved  approximately,  by  a  single  sweep  of  the  smoother,  to  provide  AilP,  Au^;  d)  Aip,  Aw  are 
evaluated  as 

(Aip.  Aw)  =  (Aip.  A<J H)  (6) 

|| 

where  IH  is  the  standard  bilinear  interpolation  operator  from  the  current  grid  H  to  the  finest  one  h;  e) 
the  (finest-grld)  solution  is  finally  updated  as 

{ip,  <*>)“♦■  {ip,  w)B  +  {Aip,  Au)  (7) 

It  is  worth  noticing  that  the  proposed  approach,  as  described  above,  is  applicable  to  any  set  of  time- 
discretized  systems  of  partial  differential  equations,  is  extremely  simple  and  does  not  require  any  “logical 
choices"  or  “free  parameters"  to  be  tuned,  insofar  as  a  single  relaxation  sweep  is  performed  at  every  grid 
level,  nor  any  artificial  dissipation  in  addition  to  that  associated  with  the  upwinding  of  the  coefficients  in 
the  LHS  of  eqns  (3,  4),  which  identically  vanishes  at  convergence.  Furthermore,  by  solving  always  for  the 
deltas  on  every  grid,  simple  homogeneous  boundary  conditions  can  be  applied  at  every  but  the  finest 
mesh,  where,  of  course,  the  appropriate  boundary  conditions  need  to  be  Imposed.  This  is  a  very  desirable 


properly  of  the  method,  especially  for  the  present  case  of  the  vorticlty  stream  function  equations,  for 
which  the  vorticlty  at  the  wall  is  obtained  by  a  Neumann  condition  for  the  stream  function,  so  that  a 
complicated  interpolation  procedure  would  be  required  if  the  vorticlty  was  to  be  computed  also  on  the 
coarser  meshes.  Moreover,  since  the  finest- grid  solution  is  immediately  updated  after  every  relaxation 
sweep,  the  correct  Neumann  boundary  condition  can  be  applied  not  only  at  the  finest- grid  relaxation 
sweeps.  Implicitly,  but  also  immediately  after  every  updating  of  the  solution,  as  an  explixit  correction.  In 
this  way  the  expected  deterioration  of  the  convergence  rate  of  the  method  in  the  presence  of  Neumann 
boundary  conditions  (an  unavoidable  trade-off  of  the  simplicity  of  the  approach)  can  be  reduced 
significantly.  Of  course,  the  extreme  simplicity  of  the  method  has  its  drawbacks.  In  particular,  the  work 
per  iteration  of  the  proposed  approach  is  slightly  greater  than  that  of  other  multigrid  schemes  and  its 
convergence  rate  is  rather  slow  at  the  beginning  of  the  process,  when  the  deltas  computed  on  the  coarser 
meshes  contain  high  frequency  errors,  which  are  thus  fed  back  into  the  sought  finest-grid  solution. 
However,  the  "pros*  are  believed  to  outperform  the  "cons",  as  it  will  be  shown  in  the  next  section. 

Results 

The  proposed  methodology  has  been  applied  to  solve  three  different  problems. 

The  Laplace  equation  in  a  unit  square  was  considered  at  first,  to  verify  the  validity  of  the  proposed 
approach  for  the  case  of  a  simple  elliptic  linear  equation  with  Dirlchlet  boundary  conditions  such  that  the 
exact  solution  of  the  differential  problem  was  given  as  exp(ffy/2)  sin(irx/2).  A  standard  line  Gauss-Seidel 
method  was  used  as  smoother  for  the  present  application,  together  with  the  collection  operator  (for  the 
residual)  called  "9-point  restriction"  by  Wesseling  (see  (17]).  Convergence  histories  or  the  method  are 
given  In  Figures  1,  where  the  (flnest-mesh)  residual  is  plotted  versus  the  number  of  Iterations.  In  Figure 
la  the  erfect  of  the  number  of  grid  levels,  the  finest  grid  being  always  a  65x65  uniform  one,  is  considered: 
the  convergence  rate  of  the  method  appears  to  Increase  significantly  with  the  number  of  grid  levels  used 
in  the  computation  and  is  similar  to  that  of  well  established  multigrid  methods.  Figure  lb  shows  Instead 
the  Influence  of  the  site  of  the  finest  grid,  when  using  always  5  grid-levels:  the  convergence  rate  is  seen  to 
be  practically  the  same  in  all  cases,  as  it  should.  The  above  results  are  considered  a  sufficient  evidence  of 
the  potential  validity  of  the  proposed  approach  and  of  the  correctness  of  its  Implementation. 

The  vorticity-stream  function  Navier-Stokes  equations,  applied  to  the  classical  driven  cavity  flow  at  Re  = 
1000,  were  then  considered  as  a  much  more  severe  test  for  the  present  method.  The  numerical  method  of 
Ref.  6,  namely  a  line  Gauss-Seidel  method  sweeping  in  alternate  directions,  was  used  as  the  smoother  in 
the  present  multigrid  strategy,  with  a  double  sweep  being  performed  at  every  grid  level,  at  each  iteration. 
The  same  collection  operator  already  used  in  the  previous  problem  and  a  nonoptlmized  unitary  time  step 
were  used  in  this  case.  Furthermore,  after  updating  the  solution  obtained  from  every  coarse-grid 
calculation,  the  vorticlty  at  the  wall  was  corrected  by  imposing  the  Neumann  boundary  condition  for  the 
stream  function  at  the  wail  on  the  finest  grid.  Only  four  grids  were  used,  because  of  the  nonlinearity  of 
the  equations  and  the  presence  of  a  Neumann  boundary  condition  for  the  stream  function.  The 
convergence  history  of  the  method  is  given  in  Figure  2,  for  four  different  uniform  (finest)  grids.  In  all 
cases  the  residual  is  seen  to  decrease  very  slowly  at  first,  but,  after  a  reasonably  small  number  of 
iterations,  the  method  takes  on  a  multigrid-type  convergence  rate.  Such  a  behavior  is  not  surprising  and 
could  have  been  easily  anticipated.  The  equations  are  in  fact  nonlinear  and  a  sequence  of  single  relaxation 
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cycles  at  every  srid  level  cannot  be  expected  to  annihilate  efficiently  the  various  frequency  errors  of  both 
the  linear  systems  (which  are  solved  incompletely)  and  of  the  nonlinear  (Newton)  Iteration  process, 
especially  since  the  solution  is  obtained  starting  from  rest.  i.  e.,  without  any  smoothing  of  the  initial  data. 
For  completeness,  the  convergence  history  of  the  smoother  (using  a  single  65x65  uniform  grid)  Is  also 
provided  in  Figure  2.  It  is  noteworthy  that  the  asymptotic  convergence  rate  of  the  multigrid  method  1s 
characterized  by  a  reduction  of  the  residual  per  iteration  equal  to  0.83,  whereas  the  corresponding  value 
of  the  single-grid  method  is  equal  to  0.S84.  The  improvement  is  substantial  even  though  the  CPU  cost  per 
iteration  of  the  multigrid  method  is  more  than  twice  that  of  the  standard  solver.  The  results  for  the 
maximum  stream  function  in  the  cavity  and  the  vortlcity  at  the  center  of  the  moving  wall,  obtained  using 
the  finest  81x81  grid,  are  finally  given  in  Table  1,  where  the  very  accurate  results  of  Ref.  1?  are  also 
given  for  comparison. 

The  lambda-formulation  Euler  equations  were  finally  considered  in  order  to  assess  the  capability  of  the 
proposed  approach  to  deal  also  with  hyperbolic  systems.  A  model  problem,  namely  the  two-dimensional 
counterpart  of  the  subsonic  spherical  source  flow  problem  studied  In  Ref.  7  was  considered  as  a  test  case. 
Solutions  were  obtained  using  the  proposed  multigrid  strategy  and  the  Block-Explicit  method  described  in 
Ref.  7,  as  *smoothera.  A  simple  injection  operator  was  used  for  obtaining  the  collected  residuals  on  the 
mesh  H  from  the  h-mesh  residuals,  because,  as  already  found  in  Ref.  20,  for  the  case  of  the  Euier 
equations,  no  advantage  was  obtained  when  using  more  complicated  colletion  operators.  Also,  the 
smoother  being  an  explicit  one  step  method,  subject  to  a  strict  CFL  condition  for  the  time  step,  a 
different  time  step  was  used  at  every  grid  level,  equal  to  the  finest-grid  time  step  multiplied  times  the 
ratio  H/h.  Results  were  obtained  using  from  one  to  four  grid  levels,  the  finest  grid  having  25x25  and 
49x49  gridpoints,  respectively.  The  convergence  histories  are  provided  in  Figures  3a  and  3b.  It  appears 
that  the  same  number  of  iterations  (i.e..»80)  are  required  for  convergence  in  both  cases,  when  using  four 
grid  levels,  a  very  typical  property  of  "sound"  multigrid  schemes;  incidentally,  the  single-grid  scheme 
requires  more  than  300  and  more  than  600  iterations,  respectively,  to  achieve  the  same  convergence  level. 
However,  considering  the  reduced  work,  3  grid  levels  appear  to  be  the  optimal  number,  especially  for  the 
first  case;  such  a  result  is  in  agreement  with  the  findings  of  Ref.  20,  where  a  similar  approach  is  used. 

Conclusions  and  future  work 

A  novel  incremental  multigrid  strategy  has  been  proposed  to  solve  the  equations  of  fluid  dynamics,  which 
provides  considerable  efficiency  gains  for  the  case  of  elliptic,  mixed  parabolic-elliptic  and  hyperbolic 
problems.  With  respect  to  current  well  established  multigrid  schemes,  the  proposed  approach  is 
conceptually  simpler  and  does  not  require  any  parameters  or  additional  artificial  dissipation.  However,  it 
requires  more  work  per  iteration  and,  at  present,  it  has  been  shown  to  be  effective  only  for  subsonic  flow 
problems  and  uniform  grids.  Future  work  is  required  to  overcome  these  two  limitations  and  to 
demonstrate  the  applicability  of  the  method  to  the  compressible  Navier-Stokes  equations  in  conservation 
form. 
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Figure  lb:  Influence  of  finest-grid  size  on  the  convergence  rate 


DRIVEN  CAVITY  PR0BLEM  (Re  =  1000.) 


ITERATI0N  NUMBER 

A  -  33x33  FlnestGrid-4Grids  B  -  49x49  FG-4G 
C  -  65x65  FG-4G  D  -  81x81  FG-4G 
E  -  65x65  FG-1G 

Figure  2:  Convergence  history  for  the  Navier-Stokes  equations 
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Table  1.  Numerical  results  for  Re  -  1000.  Driven  cavity  problem 
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Figure  3a:  Convergence  history  for  the  lambda  equations 
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Figure  3b:  Convergence  history  for  the  lambda  equations 
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